Drivers of phenotypic divergence in a Mesoamerican highland bird

Animals derive their coloration from a variety of pigments as well as non-pigmentary structural features. One of the most widespread types of pigments are carotenoids, which are used by all invertebrate taxa and most vertebrate orders to generate red, pink, orange and yellow coloration. Despite their widespread use by diverse animal groups, animals obligately obtain carotenoid pigments from diet. Carotenoid-based coloration is therefore modulated by evolutionary and ecological processes that affect the acquisition and deposition of these pigments into tegumentary structures. The Flame-colored Tanager (Piranga bidentata) is a highland songbird in the cardinal family (Cardinalidae) that is distributed from Mexican sierras through Central America up to western Panama. While female plumage throughout its entire range is predominantly yellow, males exhibit a noticeable split in ventral plumage color, which is bright orange on the West slope and the Tres Marias Islands and blood red in Eastern Mexico and Central America. We used Multiple Regression on Matrices (MRM) to evaluate the relative contributions of geographic distance, climate and genetic distance on color divergence and body differences between geographically disjunct populations. We found that differentiation in carotenoid plumage coloration was mainly explained by rainfall differences between disjunct populations, whereas body size differences was best explained by variation in the annual mean temperature and temperature of coldest quarter. These results indicate that climate is a strong driver of phenotypic divergence in Piranga bidentata.


INTRODUCTION
The variety of colors in animal integuments is perhaps one of the most striking features of the natural world. These colorations serve a variety of signaling functions, as well as non-signaling functions such as camouflage (Land & Nilsson, 2012;Bortolotti, 2006). The enormous range of animal colors we observe is derived from the complex interplay of selective light absorption by chemical pigments, structural features that affect light scattering (for example the nanostructure of butterfly wing scales and bird feathers) and other processes such as blood irrigation of exposed fleshy areas (Shawkey & D'Alba, 2017;Shawkey, Morehouse & Vukusic, 2009). Carotenoids are used by almost every vertebrate order and all invertebrate phyla to generate yellow, orange, pink and red coloration (McGraw, 2006). However, how the evolution of carotenoid-based pigmentary coloration has resulted in the colors we see at present is unclear in many taxa. Carotenoids are unique among animal pigments in that they are obtained exclusively from diet, as animals lack the metabolic routes to synthesize them de novo (Britton, 1983). In many cases, dietary carotenoids undergo chemical modifications such as ketolation before being deposited in tissues (Lopes et al., 2016;McGraw, 2006;Hill, 2000). In many vertebrates such as birds and reptiles, carotenoid-based coloration is often considered to be used as an honest signal of individual quality (Zahavi, 1975;Olson & Owens, 1998;Steffen, Hill & Guyer, 2010;Jacot et al., 2010). Since they are obtained from diet and metabolically modified before deposition, carotenoids are thought to signal a high-quality diet, as well as high immune and metabolic function in a way that cannot easily be faked (Steffen, Hill & Guyer, 2010;Johnson & Hill, 2013).
Individuals that display more elaborate signals and brighter colors are thought to incur costs associated with these signals (Faivre et al., 2003;Martínez-Padilla et al., 2010). Since carotenoid molecules readily accept electrons, they function as antioxidants and carotenoid intake is related to body maintenance (Butler & McGraw, 2010) and immune function (McGraw & Ardia, 2003); carotenoid molecules that are locked away in metabolically inactive tissue such as feathers or scales are no longer available, and their modification and deposition diverts energetic resources away from body maintenance and direct reproductive effort (Hill, 2000). In addition, increasing the elaborateness of displays increases risk of predation, and in some cases, can increase parasite loads (Baker & Parker, 1979). It is likely that decreases in ornamentation such as the limitation of carotenoid-based signals such as having red coloration on only one part of the body or having reduced color intensity have evolved as a result of trade-offs between sexually selected signaling and survival-related natural selection (Hamilton & Zuk, 1982;Aguilera & Amat, 2007;Alonso-Alvarez et al., 2008).
Unlike melanin-based coloration, where the genetic and metabolic basis of different phenotypes is relatively well understood, we know remarkably little of the mechanistic basis of carotenoid coloration outside of model organisms (Jacot et al., 2010). This complicates the task of understanding which processes underlie divergence in coloration in species that derive their coloration from carotenoids. This difference might be caused by variation of the source of carotenoids available in the environment (Grether, 2000), genetically determined differences in the metabolic conversion and pigment deposition in integument tissues (Mundy et al., 2016), or direct selection based on environmental factors such as light availability and background color (Zink & Remsen, 1986) or abundance of featherdegrading microorganisms (Burtt Jr & Ichida, 2004). The mechanisms underlying these hypotheses are not mutually exclusive, and they are expected to overlap to some extent in most systems.
This variety of mechanisms modulating color expression nevertheless result in some well-known spatial patterns such as Gloger's rule (Gloger, 1833), which predicts that endotherms living in more humid climates tend to be darker in coloration.
A study system in which closely related populations differ in the extent or intensity of their carotenoid-based pigmentary coloration would allow us to disentangle the relative contributions of these multiple selective pressures, as well as other non-selection processes such as genetic drift. Since carotenoid metabolism and deposition appear to be regulated by a large network of genes (Price-Waldman & Stoddard, 2021), if differences in color are due to fixed genetic differences between populations, we would expect to find some correlation between differences in the expressed color and genetic distance.
Similarly, changes in body size among populations can result from interactions among ecological factors such as food availability, genetic factors and other constraints such as thermoregulation efficiency. While it is intuitive that body size can evolve as a result of direct selection pressures, examples have been observed where the evidence for intraspecific changes in body size being caused by local adaptation is weak (Seeholzer & Brumfield, 2018). Both color and body size are complex traits that are affected by a number of interacting factors, and teasing apart the relative effects of climate, genetics and geographic distance is a complex problem. Simultaneous analysis of multiple sets of explanatory variables is necessary to better understand the processes that drive phenotypic divergence. With current statistical methods, we can explicitly test whether observed divergence in a trait is related to geographic distance (Isolation by distance, IBD), genetic changes, climate, or a combination of these factors.
The Flame-colored Tanager (Piranga bidentata) inhabits deciduous forests in the highlands of both of the major mountain systems in Mexico, as well as Central America and the mountains of west Panama from around 800 m.a.s.l. to the tree line (Fig. 1A). There is also an isolated population on the Islas Tres Marias off the Pacific coast of Mexico, which occupies dry shrubland and oak forest. Currently there are four recognized subspecies, described on the basis of a disjunct geographic distribution and differences in male plumage color (Howell & Webb, 1995;Hilty, 2020). This species shows marked geographic variation, with males on the Pacific side of the continent showing orange plumage while males on the Atlantic side (including Central America) show bright red coloration (Fig. 1B).
In this study, we used regression models to explore the phenotypic variation in plumage coloration and body size in Piranga bidentata and to look for relationships between this variation and environmental factors. Since this species is non-migratory exhibits largely isolated populations on a large geographic scale, we can use it as a model system to test for isolation by distance or local adaptation. We used both simple linear regression modeling and multiple regression on matrices (MRM;Lichstein, 2007) implemented in the R package ecodist (Goslee & Urban, 2007) with color distance as a response variable to test the relative contributions of geographic distance, genetic distance, distance in climate variables and differences in tree cover on the phenotypic variation in color as to distinguish between explanatory processes. We tested three different hypotheses about the factors that contribute to color variation in Piranga bidentata. First, if variation in color reflects phylogenetic relationships between subspecies, we would expect a strong relationship between divergence in color and genetic distance. Second, if plumage color is the result of local adaptation, ecological processes or phenotypic plasticity we would expect a strong relationship between divergence in color and climatic variable divergence. Finally, if color

.continued)
Upperparts dusky orange, often washed olive, heavily streaked black on back. Wings and tail blacked, edged orangish with two broad whitish wingbars and bold white spots on tip of tertials and outer rectrices.P. b. flammea: Like bidentata, but underparts paler orange. Described on basis of disjunct geographic range. P. b sanguinolenta: Like bidentata, but has head and underparts red to orange red, wingbars and tertial spots often tinged red. P. b. citrea. Like sanguinolenta, underparts bright red. Described on basis of disjunct geographic range. (C) Feather patches from which we measured color.
divergence is the result of the fixation of neutral alleles due to drift and reduced gene flow due to geographic distance, the isolation by distance model should be supported (IBD). Under this last model, we would expect a relationship between plumage color divergence and geographic distance. For body size, if climate is the proximate driver of divergence through selection on thermal constraints, we would expect a negative relationship between body size and mean annual temperature.

Genomic data
We extracted genomic DNA from 18 frozen tissue samples collected between 1991 and 2016 (P. b. bidentata: 2♂, P. b. citrea: 2 unknown sex, P. b. flammea: 1♀6♂, P. b. sanguinolenta: 4♀3♂; Table 1, Shown in a map in Fig. S1) from the MZFC tissue collection, as well as samples donated by other institutions, using the Qiagen DNeasy (Valencia, USA) kit protocol, the EPICENTRE MasterPure kit protocol as well as standard phenol-chloroform extraction. Samples deposited at the MZFC collection were obtained under a field collection permit provided by Instituto Nacional de Ecología, SEMARNAT, Mexico (FAUT-0169). We included two Piranga roseogularis and one Cardinalis cardinalis sample to act as outgroups in phylogenetic analyses. nextRad genotyping-by-sequencing libraries were obtained from our genomic DNA samples following the procedure described in Russello et al. (2015) at the SNPsaurus laboratory in the University of Oregon. Genomic DNA was first fragmented with Nextera reagent (Illumina, Inc.), which also ligates short adapter sequences to the ends of the resulting fragments. Since some of the samples had a high amount of degraded DNA, the Nextera reaction was scaled for fragmenting 20 ng of genomic DNA but 40 ng of genomic DNA was used for input to compensate, as well as to increase fragment sizes. Fragmented DNA was selectively amplified for 27 cycles at 74 degrees, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC. The nextRAD libraries were sequenced on a HiSeq 4000 platform with one lane of single-end 150 bp reads at an average depth of 54× (University of Oregon). Raw sequence reads are available at GenBank SRA (BioProject accession PRJNA757783). We assembled the resulting reads into a SNP library using ipyrad (Eaton & Overcast, 2020) using the reference approach and mapping our reads to the Northern Cardinal (Cardinalis cardinalis) reference genome (Sin, Lu & Edwards, 2020, GenBank SRA BioProject accession PRJNA642398). We filtered reads for quality by trimming reads with a PHRED score of less than 43 from the 3 end. Since many biological processes can affect the level of genetic differentiation between two samples, the clustering threshold      at which two sequences are called as orthologous cannot be assumed by default, and must be optimized empirically for each study system (McCartney-Melstad, Gidiş & Shaffer, 2019). We selected a clustering threshold of 88% based on four metrics: (1) Pearson correlation coefficient between pairwise genetic distance and percentage of missing data, (2) cumulative variance explained by the first three Principal Components, (3) total SNPs recovered, and (4) total loci recovered (McCartney-Melstad, Gidiş & Shaffer, 2019). For the final alignment, we retained loci that were present in 75% of the samples. For analyses that require unlinked loci, we further filtered our dataset using the R packages SNPfiltR (DeRaad, 2021) and vcfR (Knaus & Grünwald, 2017), retaining only loci that were more than 1000bp away from one another.

Color data
We measured the color of seven feather patches (Crown, nape, tail throat, breast, upper belly and lower belly; shown in Fig. 1C) of 109 study skins (Table 1) held at the AMNH Ornithology Collection and Museo de Zoología de la Facultad de Ciencias (MZFC), using an OceanOptics USB 2000+ spectrophotometer paired with a pulsed xenon light source. Six of these patches are carotenoid-colored, and one (tail) is melanin-colored. The measuring probe was directed perpendicularly to the feather surface and ambient light was blocked by a drilled rubber stopper surrounding the tip of the probe. We obtained spectrographic measurements of feather reflectance in the avian visual range of 300-700 nm. Feather patch reflectance was measured relative to the reflectance of a WS-2 white standard (OceanOptics) and then processed in the R package pavo (Maia et al., 2019) to smooth the spectrographic curves and correct for electric noise. From the smoothed and noise-corrected spectrographic data of each feather patch we obtained three standard variables, hue (H3), chroma (S8) and absolute brightness (B1). We used H3 instead of the more commonly used H1 because the spectral curves of carotenoid coloration show a distinct ''shoulder'' shape, always having its peak (R max ) at 700nm, showing no variation and rendering the variable uninformative. We found that calculating the peak intensity at R mid instead of R max (H3) better reflects the observed variation in carotenoid coloration. In addition, we also obtained the red chroma (SR1) of each patch, calculated as the proportion of the area under the curve in the red region (between 580 nm and 700 nm), as a complementary measure of both hue and red spectral purity. We followed the same procedure for the tail patch, which shows melanin-based coloration, to make results comparable between feather patches. We modeled the visual response to the measured color under the noise-limited visual receptor model described by Vorobyev & Osorio (1998) as implemented in the R package pavo (Maia et al., 2019) to test if any observed differences in color translate into perceptual differences in the avian visual system. The visual model we used to simulate the quantum catch was that of the Blue Tit (Cyanistes caeruleus), as it is the phylogenetically closest to our study species from models available in the package.
The source of illumination we chose was sunlight filtered through forest canopy, as it most accurately reflects the light conditions our study species is likely to be seen under in the wild.

Morphometric data
We examined 141 study skins of Piranga bidentata, held at the American Museum of Natural History (AMNH) and the Museo de Zoología de la Facultad de Ciencias (MZFC) zoological collections, covering all four-described subspecies and the majority of the species' geographical distribution (Table 1). Since juvenile males have plumage that strongly resembles that of females, we removed from analysis any specimen with yellow plumage that lacked age data that identified it as an adult (e.g., skull ossification, presence of bursa or notes on gonad development). This also allows us to mitigate the error caused by allometric changes during different stages of individual development. We used a Mitutoyo digital caliper to take external morphological measurements of standard body parts. To reduce measurement error, we used the average of triplicate measurements of each body part for analysis. We discarded the measurements of tail length because we found it difficult to replicate and highly affected by terminal wear on feathers.

Climate data
We used the R package raster (Hijmans, 2019) to extract the 19 climate variables from the CHELSA 2.1 climatology data set (Karger et al., 2017;Karger et al., 2021) in a circular 1km buffer around the geographic coordinates of each voucher specimen. For population-level analyses, we calculated a convex hull containing all specimens belonging to that population and obtained the climate variables around the hull centroid. These climatic variables correspond roughly to temperature, precipitation and seasonality variables, and many are highly correlated. To reduce model overfitting and mitigate the risk of multicollinearity, we avoided using variables with more than |0.7| correlation in the same model as suggested by Dormann et al. (2013). To test whether there was a correlation between feather color and vegetation cover, we obtained an estimate the percentage tree cover from the Moderate Resolution Imaging Spectroradiometer (ModIS) satellite though their MOD44B.006 Terra Vegetation Continuous Fields Yearly 250 m dataset (DiMiceli et al., 2015), accessed through the Google Earth Engine portal (Gorelick et al., 2017) and operated by NASA. We calculated the median vegetation cover from 19 yearly images spanning the period from 2000 to 2018. Finally, we calculated pairwise geographic distance between each collection point and between population centroids using the great circle distance algorithm implemented in the R package geodist (Padgham & Sunmer, 2019).

Analysis
To first assess if genetic and color variation are phylogenetically structured (color groups corresponding to described subspecies), we constructed a phylogenetic tree with our SNPs data. We used the SVDQuartets (Chifman & Kubatko, 2014) algorithm implemented in PAUP* (Swofford, 2002) to obtain a coalescent tree topology. We generated 100 bootstrap replicates to assess branch support in the resulting consensus phylogenetic reconstruction, which we present in Fig. 2. We then assessed visually if members of the described subspecies clustered together into monophyletic groups. We also assessed these groups using the Bayesian clustering algorithm implemented in the software ADMIXTURE (Alexander, Novembre & Lange, 2009), using a K determined by minimizing cross-validation error. This allowed us to estimate the ancestry of individuals within our sample. As the study species is obviously dichromatic, we performed separate analysis on colorimetric data by sex. We explored variation in color by performing one-way analysis of variance (ANOVA) on the calculated color variables (B1, S8, H3 and S1R) with the subspecies as a factor. The results of this analysis are in the Tables S1 and S2. To explore if intraspecific variation is consistent with Gloger's Rule we constructed linear models with plumage brightness (B1) as a response variable and mean annual precipitation (bio12) and mean annual temperature (bio1) as predictors. Ventral coloration is under stronger sexual selection compared to dorsal coloration in many bird taxa (Shultz & Burns, 2017;Marcondes & Brumfield, 2019). Since ventral color is relatively uniform patterned in this species, we took the average of the four ventral patches as representative of overall carotenoid color in subsequent analyses.
To further assess if genetic clusters match color groups, we performed linear discriminant analyses on ventral coloration variables on females (using plumage brightness, hue and chroma; B1, H3 and S8) and males (using plumage brightness, hue, overall chroma and red chroma; B1, H3, S8 and S1R) separately.
We followed a similar procedure for morphometric data, using one-way ANOVA to test for differences in each measured variable between subspecies (Table S3). Since one morphological variable was one order of magnitude larger than the rest, we rescaled all morphological data by log -converting it. To test if there were differences in body size among the described subspecies, we carried out an analysis of variance (ANOVA) with subspecies as a factor for each of the variables we measured. The results of these tests are shown in Table S3. To condense this morphological variation, we carried out a principal component analysis (PCA) on our morphological data matrix and obtained a set of principal components, the first of which explains 79.4% of the total variation, and is positively correlated almost entirely with wing chord. The second principal component explains 12.76% of variance and is positively correlated with the rest of the morphological variables. To test if this species follows the pattern known as Bergmann's Rule, we constructed a linear regression model with the first PCA score value as a response variable and mean annual temperature (bio1) as the predictor.
To distinguish between the effect of genetic distance, geographic distance, and environmental variables, we used multiple regression on matrices (MRM) as implemented in the R package ecodist (Goslee & Urban, 2007), which tests for relationships between two or more distance matrices using matrix permutation. We also tested the relative effects of geographic distance and rainfall divergence on allele frequencies using a Bayesian MCMC framework implemented in the R package BEDASSLE (Bradburd, Ralph & Coop, 2013). We ran BEDASSLE for 1,000,000 generations using the beta binomial model, which allows populations to diverge from model expectations due to overdispersion. We used a matrix of Euclidean distances on rainfall, calculated from the CHELSA climatology data set (Karger et al., 2021). To test for the effects of selection, we used Bayescan 2.1 (Foll & Gaggiotti, 2008) to identify loci under selection in the color morphs. These tests allowed us to test the likelihood of the null model of isolation by distance (IBD), divergence due to drift, and local adaptation. Since we lack a one-to-one match between genomic sequences and color data, we carried out this analysis using the average spectrum per subspecies and an Fst distance matrix between the four subspecies.

Genomics
From our data set of 84,739 SNP loci we obtained a coalescent tree topology (shown in Fig. 2A) which grouped all of the Piranga bidentata samples into a single well-supported clade, supporting the monophyly of Piranga bidentata. The two western groups with orange male plumage (bidentata and flammea) cluster into one subclade with middling support. Samples from populations with red male plumage cluster together, although the topology doesn't reflect the geographic structure of the samples, and we did not recover the subspecies P. b. sanguinolenta as monophyletic. Despite our limited sampling, we found some evidence for genetic structure in this species.
Our ADMIXTURE results (K = 2, Fig. 2B) show that most of our samples cluster unambiguously into two groups. These groups correspond to the orange color morph (subspecies bidentata and flammea) and the red color morph (subspecies citrea and sanguinolenta). Our two samples from the bidentata subspecies appear in the analysis as admixed between these two groups.
We performed two Bayescan runs of 100,000 generations each on a subset of 8,440 unlinked loci, one with prior neutral model odds of 100 and the other with odds of 1,000. In both runs we failed to detect any loci as statistically significant outliers in Fst, suggesting the markers in our data set are selectively neutral. However, this method risks low statistical power with small sample sizes, so these results should be interpreted carefully. A more extensive sampling is needed to fully assess selection in this taxon.
Our BEDASSLE run resulted in an aE/aD posterior distribution with a mean of 15,091 which suggests that the effect of rainfall on allele frequencies is relatively small compared to the effect of geographic distance (A difference of one mm of annual rainfall having a comparable effect as about 15km of lateral movement).

Color
The ANOVA of the color data (Table S1) showed a break in males' carotenoid coloration for all four of the variables we computed (H3, S8, SR1 and B1), with males of the orange group (subspecies bidentata and flammea) having brighter, more orange and less spectrally saturated plumage color than males from the red group (subspecies sanguinolenta and citrea).
Our LDA result shows that for males (Fig. 3A), ventral color variables B1, H3, S8 and S1R distinguish two color groups congruent with the genetic clusters we detected. However, these predictors can't distinguish between the described subspecies that form each of the groups. We found a similar pattern for females (Fig. 3B), although the separation is much less clear.
Under the noise-constrained visual receptor modeling (Vorobyev & Osorio, 1998), the average difference in color that we measured between the red and orange groups is higher than the perceptual threshold (i.e., larger than 1 JND, for just-noticeable difference Fig. 4A), which indicates that the colors are sufficiently different that birds are capable of discriminating between them. Notably, there was also a slight but noticeable and statistically significant difference between the West mainland and the Tres Marias orange birds, with the island birds having brighter orange plumage. The pairwise differences in ventral color between subspecies of the same color group (bidentata-flammea and sanguinolenta-citrea) are much smaller. . Linear discriminants were calculated from the variables B1, H3 and S8 for both, with the addition of S1R for males. For males, we can observe discrimination between the orange (bidentata and flammea) and red (citrea and sanguinolenta) groups, but not between the subspecies comprising each group. For females, we can see a similar pattern, but less strongly differentiated.
Full-size DOI: 10.7717/peerj.12901/ fig-3 Females were generally consistent in their ventral carotenoid coloration, showing no significant variation in either hue or chroma among subspecies (Table S2). We found a small but statistically significant variation in plumage brightness among subspecies, which follows the same overall pattern as in the males. Although the pairwise differences in overall carotenoid color between subspecies are larger than one 1 JND, the threshold falls within the confidence interval (Fig. 4B).
For melanin coloration, we found a similar pattern and magnitude of divergence between subspecies in both males (Fig. 4C) and females (Fig. 4D).
For males, we found a significant negative relationship between plumage brightness (B1) and annual precipitation (bio12), for both carotenoid (ventral, F 1-73 = 26.6, adj. R 2 = 0.26, p > 0.001, Fig. 5) and melanin (tail, F 1-73 = 13.62, adj. R 2 = 0.15, p < 0.001) feather patches. We found a similar negative trend in female color, although the coefficients of determination were lower (ventral, F 1-48 = 12.62, adj.R 2 = 0.19, p < 0.001; tail, F 1-47 = 6.34, adj. R 2 = 0.10, p = 0.015). These results show that this species follows the ecogeographical pattern known as Gloger's Rule, which predicts that among ectothermic species, individuals found in habitats with higher humidity will tend to have darker coloration. We found a similar negative relationship between ventral plumage brightness and isothermality bio3 (males, F 1-73 = 27.42, adj. R 2 = 0.26, p > 0.001; females, F 1-48 = 8.76, adj. R 2 = 0.14, p = 0.005), as well as a statistically significant but weak negative relation between tail brightness and isothermality (F 1-73 = 7.22, adj. R 2 = 0.08, p = 0.009) in males. To further test whether there was support for the hypothesis that Gloger's Rule is at least partially due to selection of darker individuals in habitats with more dense vegetation, we tested for a relationship between plumage brightness and remotely sensed vegetation coverage calculated from ModIS measurements. We found a significant relationship between plumage brightness (F 1-60 = 16.55, adj. R 2 = 0.2, p < 0.001) and vegetation cover, but not for red chroma (F 1-60 = 1.37, adj. R 2 < 0.001, p = 0.246) and vegetation coverage, as calculated from ModIS measurements. The results of MRM analysis showed that for males there is statistically significant positive relationship between differences in breast feather color and geographic distance (F = 73.26, R 2 = 0.037, p = 0.001), and no significant relationship between breast color differences and genetic distance (F = 1.97, R 2 = 0.33, p = 0.38) or difference in vegetation cover (F = 9.26, R 2 = 0.004, p = 0.26). However, we found a significant positive relationship between differences in mean annual rainfall and differences in breast color (F = 439.229, R 2 = 0.19, p = 0.001). For females, we found no Figure 5 Negative relationship between ventral plumage brightness and annual precipitation for males. This relationship is consistent with the ecogeographic pattern known as Gloger's rule, which predicts that organisms living in more humid climates will tend to be darker. Full-size DOI: 10.7717/peerj.12901/fig-5 relationship between color distance and geographic distance (F = 16.064, R 2 = 0.015, p = 0.128) or differences in rainfall (F = 1.92, R 2 = 0.002, p = 0.544). The MRM model with male breast color distance as response variable and Fst, geographic distance and rainfall distance matrices as joint predictors shows a not statistically significant but very large proportion of color variance is explained by these three factors (F = 26.39, R 2 = 0.975, p = 0.07). We obtained a similar result with female breast color distance, although the trend is much weaker (F = 0.50, R 2 = 0.43, p = 0.79)

Morphometrics
We found no significant sexual dimorphism in this group for five of the six morphological variables we measured. The only measurement for which sexual dimorphism was statistically significant was the wing chord. Of the six measurements, we found significant variation between subspecies in 5 of them (Table S3). The first two Principal Components resulting from our PCA explained 73.48% of the observed variation. Principal Component 1 was strongly positively correlated with wing chord, and Principal Component 2 was correlated with the other morphological variables we measured, and thus serves as a proxy for overall body size. We found that PC1 scores separate the subspecies bidentata from the rest, while PC2 separates the Mexican subspecies (bidentata and sanguinolenta) from the remaining two. However, all subspecies showed a high degree of overlap in morphometrics (See Supplemental Information). MRM results showed a small but significant relationship between our proxy for body size (PC2) and geographic distance (F = 202.57, R 2 = 0.02, p = 0.001) but not between PC2 and mean annual temperature (F = 0.3, R 2 = 0.0003, p = 0.84).

DISCUSSION
Our results show that there is significant variation in male feather coloration and body size among populations of Piranga bidentata, as well as evidence of genetic structure. This variation in color covaries in different amounts with genetic structure, geographic distance and climatic variables. These results suggest support for the hypothesis that the observed variation in color is partially the result of local adaptation, phenotypic plasticity or a combination of the two. The most striking of these differences is the change in male carotenoid coloration from orange in the Western Mexico subspecies (bidentata and flammea) to red in the Eastern Mexico and Central American subspecies (citrea and sanguinolenta), which we found to be related to mean annual rainfall. Overall plumage brightness also differs among groups, and this variation is related to mean annual rainfall, but not to vegetation coverage calculated from ModIS measurements. These results are consistent with the well-known ecogeographical pattern known as Gloger's rule (Gloger, 1833). For body size, we found no relationship between a measure of body size and mean annual temperature, but there is a significant relationship between body size and geographic distance, suggesting that the mechanism at work is isolation by distance. A genomic scan for loci associated with body size could clarify the role of genetic differentiation on variation of body size.

Plumage color
Plumage carotenoids are obtained exclusively from the diet, but the relationship between carotenoid ingestion and plumage color expression is not necessarily direct. The deposition of red and orange carotenoids (e.g., canthaxanthin and astaxanthin) requires metabolic conversion through ketolation of dietary yellow carotenoids (e.g., zeaxanthin). Although the metabolic pathway is not fully understood, evidence from several studies suggests that CYP2J19, a member of the P450 cytochrome gene family, is implicated in the ketolation reaction of dietary carotenoids (Mundy et al., 2016;Lopes et al., 2016;Emerling, 2018;Sin, Lu & Edwards, 2020). The enzymatic activity of CYP2J19 and other cytochrome proteins appears to be highly correlated with overall mitochondrial activity (Hill et al., 2019;Cantarero et al., 2020), which potentially maintains the honesty of modified carotenoid ornamentation as a signal of individual quality. While Koch et al. (2019) found that the evidence linking carotenoid deposition to several measures of individual performance is weak, other studies show evidence that carotenoid pigmentation is correlated with testosterone levels (Khalil et al., 2020).
While the exact chemical makeup of the feather carotenoids in this species has not yet been studied, closely related species of the genus Piranga derive their yellow feather coloration from dietary zeaxanthin and their red and orange coloration from a combination of several 4-keto-carotenoids such as canthaxanthin and astaxanthin (Hudon, 1991;Lopes et al., 2016). This seems to be largely conserved in the genus, with a few exceptions (Hudon, 1991). Since the difference in color that we observed appears to be continuously variable, it is likely that observed differences in color result from changes in the relative proportions of several carotenoids being deposited in feathers. This ratio is likely associated to some degree with the activity of CYP2J19-like gene products, which perform these ketolation reactions.
Another potential source of variation in coloration is the effect of differential sexual selection between populations. Our pairwise comparisons (Fig. 3) show similar patterns and magnitude of variation between subspecies for female carotenoid (3B), male melanin (3C) and female melanin (3D), but much larger differences in male ventral carotenoid colors (3A). As ventral coloration has been found to be under heightened sexual selection compared to other plumage patches (Shultz & Burns, 2017), this greater variation could be the result of stronger sexual selection for ketocarotenoid coloration in the red group (sanguinolenta and citrea). A series of studies on foraging behavior and characterization of dietary carotenoids available to this species would help illuminate other possible mechanisms modulating color expression.

Gloger's rule
We found that in Piranga bidentata, darker and redder coloration in males is correlated with higher environmental humidity (measured as annual rainfall), compared to orange coloration. This is consistent with the ecogeographical pattern known as Gloger's rule, which in its original formulation predicts that for endothermic organisms, darker coloration is found in hotter and more humid environments (Gloger, 1833). This pattern is widely observed in birds and mammals, and a study by Zink & Remsen (1986). found that a majority of North American bird species studied conform to it to some degree.
One of the most widely accepted hypotheses to explain Gloger's rule is selection for darker and more cryptic plumage in environments with dense vegetation and thus lower environmental light levels (McNaught & Owens, 2002). Alternative explanations involving other, non-mutually exclusive processes have also been proposed. These include increased feather resistance to mechanical wear (Bonser, 1995), increased resistance to bacterial and fungal degradation of feathers (Burtt Jr & Ichida, 2004;Gunderson et al., 2008), and pleiotropic effects from selection on other traits which increase with humidity such as immune responsiveness and overall mitochondrial activity (Johnson & Hill 2013;Hill et al., 2019;Cantarero et al., 2020). In general, these proposed explanations suggest that Gloger's rule results from an increase in pigment load in feathers, stemming from selection on other traits. We observed that plumage brightness of both melanin and carotenoid-colored patches is negatively correlated with rainfall (though not with a measure of vegetation cover) as predicted by Gloger's Rule, similar to what has been observed in Australian birds (Delhey, 2018) and in other American birds (Ramírez-Barrera et al., 2019). However, this contrasts with the results of Ribot et al. (2019), who found a strong association between plumage brightness and remotely sensed vegetation cover. Taken together, these results fail to support the hypothesis that darker plumage in this species is linked to a darker environment caused by denser vegetation coverage. An alternative explanation that is compatible with our results would be that the mechanism in play is selection on stronger sexual signaling and not increased crypsis. Somewhat counterintuitively, additional pigment load increases spectral purity (chroma) but decreases overall brightness, since pigments function by selective absorption of incoming light. As ambient light decreases, more pigment is necessary to achieve the same level of background contrast, but it also decreases the overall amount of light being reflected. Our data does not support the hypothesis that differences in plumage color are caused primarily by fixed genetic factors. A mutation which results in reduced function of CYP2J19-like genes or other genes involved in the carotenoid metabolism pathway would be unlikely to become fixed outside of founder effect scenarios, as these mutations appear to occur in very low frequencies in the wild and would be rather strongly selected against. We found that the color of males from the Tres Marias islands is significantly different from mainland males, and that this difference is strong enough to be distinguishable by avian visual systems. The color of Tres Marias males is achromatically brighter and less saturated compared to their mainland counterparts. This is consistent with the overall pattern we found, as the climatic conditions on the islands are more xeric than on the mainland, and the dominant vegetation presents a more open canopy. However, many island species across many taxa present reduced sexual signaling when compared to their mainland sister groups, and it has been hypothesized that this is due to changes in the intensity of sexual selection compared to the mainland.
We found no significant relationship between plumage brightness or red chroma and vegetation cover, which also shows limited support for the hypothesis that increase of red coloration is due to an increase in dietary carotenoid availability linked to the environment's primary productivity as was the case in the work of Jacot et al. (2010). The strong positive correlation between feather color and rainfall suggests that the mechanism modulating carotenoid deposition is related to humidity, such as selection acting on immune response or mitochondrial activity. A direct measurement of mitochondrial activity or quantitative transcription of genes associated with immune response would allow us to directly test this hypothesis.
The lack of relationship between color variation and geographical distance indicates that the mechanism underlying this variation is not isolation by distance. This is in contrast to the findings of Núñez Zapata et al. (2018) and by Ramírez- Barrera et al. (2019), who used a similar approach and found that for Turdus assimilis and Habia rubica respectively, Isolation by Distance is the model that best explains color variation. This, coupled with the lack of genetic structure that correlates to geography suggests that there may be more gene flow than previously thought and/or that carotenoid-based coloration in this species is fairly plastic.
Coloration is a complex and likely highly modular phenotypic trait. While carotenoidbased coloration is undoubtedly under some degree of fixed genetic regulation, both experimental (e.g., McGraw & Hill, 2000;Jacot et al., 2010) and comparative studies such as Ramírez- Barrera et al. (2019) and this one show that it is relatively plastic and subject to a large degree of environmental modulation. Our results show that a single gross-scale environmental measurement like mean annual rainfall can have a disproportionate effect on the coloration of certain taxa. However, the exact mechanisms that mediate the observed color differences have yet to be tested, and experimental studies that explicitly test these hypotheses are needed.

Body size
Our results on body size are not conclusive; we found a statistically significant relationship between body size and geographic distance, suggesting that the mechanism driving body size variation is isolation by distance. This is similar to what was found in Cranioleuca antisiensis, where the mechanism underlying a dramatic change in body size between populations was found to be isolation by distance (Seeholzer & Brumfield, 2018). Similar results were also found by Ramírez- Barrera et al. (2019) in Habia rubica. These studies make apparent that variation in body size in these neotropical birds is mostly driven by geographic isolation. However, although the relationship we found was statistically significant, the effect size was very small and we would be cautious about extracting strong conclusions from it. A genomic scan of alleles linked to variation in body size would further clarify the role of genetic variation on phenotypic variation this species.
Taken together, our results show that coloration and body size evolution are seemingly decoupled in this species. In conclusion, we found that climate is an important driver of phenotypic divergence in this species. An experimental study manipulating availability of dietary carotenoids would help illuminate the proximate mechanisms underlying this divergence. Although our phylogenetic analysis shows that the two samples from Central America (ssp. citrea) separate into a well-supported clade, a more extensive sampling is required to clarify their taxonomic status.
Finally, although we find evidence of genetic structure in this taxon, we suggest that a broader sampling and more in-depth genomic analysis is necessary to fully assess the relationship between genomes, environment and phenotypes.